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ABSTRACT 



The periodogram is a popular tool that tests whether a signal consists only of noise or if it also includes other compo- 
nents. The main issue of this method is to define a critical detection threshold that allows identification of a component 
other than noise, when a peak in the periodogram exceeds it. In the case of signals sampled on a regular time grid, deter- 
mination of such a threshold is relatively simple. When the sampling is uneven, however, things are more complicated. 
The most popular solution in this case is to use the Lomb-Scargle periodogram, but this method can be used only when 
the noise is the realization of a zero-mean, white (i.e. flat-spectrum) random process. In this paper, we present a general 
formalism based on matrix algebra, which permits analysis of the statistical properties of a periodogram independently 
of the characteristics of noise (e.g. colored and/or non-stationary), as well as the characteristics of sampling. 

Key words. Methods: data analysis - Methods: statistical 



1. Introduction 

Spectral analysis is a popular tool for testing whether a 
given experimental time series {x{tQ),x{ti), . . . ,x{tM-i)} 
contains only noise, i.e. x(tj) — n{tj), or whether some 
other component s{t) is present, i.e. x{t) = s{t) + n{t). The 
classic approach is to fit the time series with the model 
function 

N-l 

x{tj) = COS {2n fktj) + bk sin (27r/feij) + n{tj), (1) 

k=0 

J = 0, 1, . . . , Af — 1. If, for example, a periodic component 
s(t) is present with a frequency /; close to one in the set 
{fk}, then the periodogram {pk}, 



Pk = al + bl, 



fc = 0, 1, 



,iV-l, 



(2) 



will show a prominent peak close to k = I. If s(t) is semi- 
periodic or even non-periodic, the situation is more com- 
plicated since more peaks are expected. The main problem 
with the use of this technique is the definition of a detection 
threshold that fixes the contribution of noise in such a way 
that, when a peak exceeds it, the presence of a component 
s{t) can be claimed. In the case of signals sampled on a 
regular time grid, the determination of such a threshold is 
a relatively simple procedure, but this is not the case when 
the condition of regularity does not apply. I n this respect , 
several solutions have been propo s ed (see Lomb I 1975 
Mello I 119811: IScargle I [198^ iGiUiland fc Baliunas I 



Il987t iReegen |[2007t IZechmeister fc Kiirster Il2009l and ref- 
erences therein) that, however, work only under rather re- 
strictive conditions (e.g. white and/or stationary noise) and 
are difficult to extend to more general situations. 

In this paper, a general formalism is presented that al- 
lows analysis of the statistical properties of periodograms 
independently of the specific characteristics of the noise and 
the sampling of the signal. In Sec. [5] the formalism is pre- 
sented for the case of even sampling and its extension to 
arbitrary sampling in Sec. [31 The usefulness of the proposed 
formalism is illustrated in Sec. |31 where the case of white 
noise with a mean different from zero and that of colored 
noise is calculated. In Sec. [S] the relationship between the 
periodogram and the least-squares method is considered. 
Finally, on the basis of simulated signals and a real time 
series, we discuss in Sec. [S] whether the use of algorithms 
specifically developed for computing the periodogram of 
unevenly-sampled time series is really advantageous. 



2. Periodogram analysis in the case of even 
sampling 

If a continuous signal x{t) is sampled on a set of N eq- 



uispaced time instants to,ti, 
j = 0, 1, . . . , — 1 is obtained, 
form (DFT) is given by 



Xk 



erraz- 



N-l 

E 

J=0 



.jtjsf-i, a time series x 
Its discrete Fourier trans- 



fc = 0, l,...,iV- 1, (3) 



2 
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with i = \/— 1. The sequence {xj} can be recovered from 
{xk} by means of 



1 

k=0 



Xke 



2-iTkj/N 



j = 0,l,...,iV-l. (4) 



The set {k/N}^^^^ provides the so-called Fourier frequen- 
cies. Implicit in the use of DFT is the assumption that 
{xj} is a periodic sequence with period A^Ai where Ai = 

In matrix notation, Eqs. ([3]) and (|4]) can be written in 
the form 

x = Fx (5) 

and 

X = F*x. (6) 

Here, x and x are column arrays that contain, respec- 
tively, the sequences {xj} and {xk}, and F is the so-called 
^^Fourier matrij^' , which is an A'^ x iV square, symmetric 
matrix whose (A:, j)-entry[3 is given by 



F, 



kj 



-i2-nkj/N 



(7) 



The superscript "*" denotes the complex conjugate trans- 
pose. Matrix F is unitary, i.e., 

FF* = F*F = J, (8) 

with I the identity matrix. Another useful property is that 



FF = F'^ F = FF'^ 



H 



with 



H 



By means of Eq. 

F^Fx - 



/ 1 ••• 0\ 
••• 10 
••• 10 



1 ••• 
VO 1 ••• 0/ 

it can be shown that 



(9) 



(10) 



F'^Fx 



F^Fl 



0, 



(11) 



where Fsr = 5R[F] and Fx = T[F] and where and I[] 
are the real and the imaginary parts of a complex quantity. 
Indeed, 



FF = Fk-Fsr + FxFx + i2F^Fx = H, 



(12) 



and i? is a real matrix, hence Eq. (TTTT) has to hold. From 
Eq. ([3]) it is also possible to see that xq is a real quantity. 
In the case that N is an even number, the same holds for 
a;Ar^+i where = \N/2] and [2;] represents the smallest 
integer greater than z. Finally, dealing with x, only half 
of this array can be considered, since xjq-k — ^*ki k — 
1, 2, . . . , A^l ~'^27V-f ' where 51 — \ \ii — j and zero otherwise. 
With this notation, the periodogram of x is defined as 



Pfe = 2(5R[x,]2+I[x,]2) = 2|£,| 



,k\ , fc = 0,l,...,A^t-l' 
(13) 



^ In the following, the element in the nth row and mth column 
of an Ai' X M matrix A will be indicated with Amn or alterna- 
tively with (A)m„, n = 0, 1, . . . A - 1, m = 0, 1, . . . , M - 1. 



where 



is the fcth entry of the array S — 



[xo, xi, . . . , Jcat^-i] and |.| denotes the Euclidean norm. 
By means of 



(14) 



a column array obtained by the column concatenation of 
Sjj = 3fi[S] and = Eq. can be rewritten in the 
form 



Pfc = 2(Sf + zL+fc), fc = 0,l,. 



,iVt 



1. 



(15) 



In Sect. [S] it is shown that this periodogram is equivalent to 
that obtainable by means the least-squares fit of model ([ij 
with M = A^, tj = j and fk = k/N. 

An important point to stress is that, if x is the realiza- 
tion of a (not necessarily Gaussian) random process, then 
each Xj. is given by the sum of N random variables. This is 
because of the linearity of the Fourier operator F . Thanks 
to the central limit theorem^ therefore, the entries of x can 
be expected to be Gaussian random quantities. As a conse- 
quence, the entries of z can also be expected to be Gaussian 
random quantities with covariance matrix = E[££"^] 
given by 



FxCxF^ 



Fl 771 jpl 
sc. £_3i'~^x£_X 

FjtCxFt 



(16) 



Here, E[.] denotes the expectation operator^ Cx ~ E[xx'^] 
is the covariance matrix of x, and £ the matrix obtained by 
the first rows of the Fourier matrix F. From Eqs. (1161) 
and pT|) . it is easy to deduce that, if x is the realization 
of a standard white- noise process, i.e., Cx = I, then Cs 
is a diagonal matrix with (C£)ii = 1, {C£)NfN^ — and 
{Cz)kk = 0.5. In other words, the entries of z are mutually 
uncorrelated. In turn, this means that 5J[xj,] is uncorrelated 
with If a; is a colored (not necessarily stationary) 

noise process, i.e. Cx is not a diagonal matrix, then this 
holds also for C^. However, from x it is possible to ob- 
tain an array y containing mutually uncorrelated entries 
by means of the transformation 



y 



1/2 

The matrix Cx can be computed via 



c: 



-1/2 



-1/2 



u 



(17) 



(18) 



with U the orthogonal matrix whose columns contain the 
eigenvectors of Cy and S a diagonal matrix containing the 
corresponding eigenvalues {A;}^q^. This decomposition is 
particularly simple and computationally efficient if Cx is a 
circulant matrix since it can be diagonalized according to 



F*diag[Fc]F, 



(19) 



where "diag[q]" denotes a diagonal matrix whose diagonal 
entries are given by the array q and c is the first column 
of Cx. Because of this, the decomposition (IT51) can be di- 
rectly computed with U = F and S = diag[c]; hence, 
y = F*'E~^^^F. This means that the whitening operation 



^ From now on, if r is an A x 1 column array, then r is a 
column array that contains the first Af = [A/2] entries of r, 
i.e. r = [ro, ri, . . . , rN.f-i]'^ . Similarly, if A is an A x M matrix, 
then A is a matrix that contains the first Af rows of A. 
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can be performed in the harmonic domain through the fol- 
lowing procedure: a) computation of x, i.e. the DFT of a;; b) 

^ ^ 1/2 

computation of {yk} = {xk/^i/ }; and c) the inverse DFT 
of y. A potential difficulty in using transformation ([T7| is 
that Cx is required to be of full rank. Such a condition can 
be expected to be satisfied in most practical applications. 
Some problems can arise if the time step of the sampling 
is much shorter than the decorrelation time of nE|. Indeed, 
some columns (and therefore some rows) of C could almost 
be identical i.e., this matrix could become numerically ill- 
conditioned. In this case, the most natural solution consists 
of averaging the data that are close in time. 

When the periodogram is used to test whether x — n 
vs. X = s + n, a. threshold Lp^ has to be defined such that, 
with a prefixed probability, a peak in pk that exceeds Lp^, 
can be expected to not arise because of the noise. This re- 
quires knowledge of the statistical properties of p under the 
hypothesis that x — n. The simplest situation is when n 
is the realization of a standard white-noise process. In fact, 
since the entries of z are uncorrelated random Gaussian 
quantities, from Eq. (jl5p it can be derived that the en- 
tries of p are (asymptotically) independent quantities dis- 
tributed according to a X2 distributiorQ. As a consequence, 
independent of the frequency fc, a threshold Lpa can be de- 
termined that corresponds to the level that a peak due to 
the noise would exceed with a prefixed probability a when a 
number Nf of {statistically independent) frequencies are in- 
spected. More specifically, Lpa is the highest value for which 
1 - [1 - exp {-LpJ]'^f < a fScargle 1982), in formula 



Lpa = sup {1 - [1 - exp {-LpJ]'^f < a} 



(20) 



If the entire periodogram is inspected, then Nf = N^. 
Commonly, Lpa is called the level of false alarm. 

Threshold (|20p is not applicable when the noise is col- 
ored. However, the requirement to fix a different level for 
each frequency can be avoided if the original signal x is 
transformed into 

y = c;,'^'x. (21) 

Indeed, under the hypothesis that x = n, the entries of y 
are uncorrelated and unit-variance random quantities. As 
seen above, this operation should not be difficult. Some 
problems emerge if the decomposition (llSp is carried out 
by means of the efficient DFT approach described above 
(a necessary approach in the case of very long sequences of 
data). This is because, in forming C„, it is necessary to take 
the periodicity of the sequence x forced by the DFT into 
account. In other words, it is necessary to impose a spu- 
rious correlation among the first and the last entries in n. 
For instance, for a stationary noise process, Cn is a Toeplitz 
matrix, but it has to be approximated with a circulant one. 
When dealing with sampled signals, this is an unavoidable 
problem that no technique can completely solve. A classical 
solution to relieving this situation is the windowing method^ 
i.e., the substitution of Xj in Eq. with rijXj, where {rjj} 
is some prefixed discrete function (window) that makes the 
signal gently reduce to zero at the extreme of th e sam- 
pfing interval (e.g. see lOpoenheim fc Shafer I Il989t ). This 



^ The decorrelation time of a random signal n{t) is the time 
interval At such that two values n{t\) and n(t2), with t2 ~t\ = 
At, can be considered as uncorrelated. 

* xi denotes the chi-square distribution with two degrees of 
freedom. 



method can be expected to work satisfactorily only when 
the decorrelation time of noise is shorter than the length of 
the signal. 

3. Periodogram analysis in the case of uneven 
sampling 

If a signal x{t) is sampled on an uneven set of time instants, 
some problems emerge: it is no longer possible to define 
a set of " natural frequencies such as those obtained by 
the Fourier transform. In turn, this implies some ambigui- 
ties in the definition of the Nyquist frequency that, loosely 
speaking, corresponds to the highest frequencie s that con- 
tain informati on o n the signal of interest (e.g. see lVio et al. I 
I2OOOI : lKoen ll2006l) . As a consequence, Eq. Q has to be 
modified. In the following, with no loss of generality, it is 
assumed that x{t) is sampled at M arbitrary time instants 
to, ti, . . . , tM~i with to — 0, tM-i = M —I and the remain- 
ing tj arbitrarily distributed within this interval. Moreover, 
a set of N frequencies A: = 0,1,...,A^— lis considered with 
N ^ M. Such a set corresponds to the frequencies that are 
typically inspected when looking for a periodicity. However, 
others can be chosen. With these conditions, a transforma- 
tion corresponding to the one given by Eq. ([5]) is 



X = Tx, 



where 



^kj - 



-i2-iTktj /N 



M 



(22) 



(23) 



ij = tj/ Amt, Ajnt — 7min [{ij+i — tj}] and 7 is a real pos- 
itive number. Because of the normalization by A,„t, time 
tj is expressed in units of the shortest sampling time inter- 
val. Apart from the substitution of the Fourier matrix F 
with J^, the uneven sampling of signals does not modify the 
formalism introduced in the previous section. In particular, 
the covariance matrix Cg defined in Eq. (jl6|) becomes 



(24) 



The N X AI Fourier matrix ^ does not have the prop- 
erties ([51)- (fT^ . As consequence, and also in the case that 
X is the realization of a standard white-noise process i.e. 
Ct = I, matrix Cj, 



(25) 



is not diagonal. In general, is not even diagonalizable. 
For example, this happens when the periodogram of a time 
series containing M data is computed on N frequencies 
with N > M (a, typical situation in practical applications) . 
This implies that the entries of z cannot be made mutually 
uncorrelated. Obviously, the same holds for the entries of p 
as given by Eq. ([T5|) . As a consequence, although a number 
N of frequencies are considered in S^, at most only M/2 of 
them are statistically independent^. Particularly trouble- 
some is that, for a given frequency k, ^[x_j^], and I[x_^.] (i.e. 



This is because, if > M, the rank of the N x N matrix 
is smaller than or equal to M. This implies that the array 
z has at most M degrees of freedom. Since each entry of p is 
given by the sum of two entries of z, then a periodogram has at 
most M/2 degrees of freedom. 



4 
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Zfc and zjvt+fc) are also correlated. This makes it difficult 
to fix the statistical characteristics of pt- In this respect, 
two choices are possible. The first consists in the deter- 
mination, for each frequency, of the PDF of pk- Actually, 
this is a rather involved approach, because Zfc and ZTv^+fc 
have variance {C£)kk and {C£)N^+k,N^+k, respectively, and 
covariance {C£)k,N^+k- Therefore, once 'zk and 'zpj^^+k are 
normalized to unit variance, each pk is given by the sum of 
two correlated random qua ntities. Although available in 
analytical form (iSi mon 1 I2006D . the resuhing PDF is rather 
complex a nd hence difficu lt to handle (for an alternative ap- 
proach, see lReegen Il2007f) . Moreover, there is the additional 
problem that Lp^, changes with k. A simpler alternative is 
the use of two uncorrelated and unit- variance random quan- 
tities, Vk and VN^+k, obtained through the transformation 

d = ^-^'^U,z,. (26) 

Here, v = [vk,VN^+k], = [zk,ZN^+k], S~ ' is a diago- 
nal matrix whose entries are given by the reciprocal of the 
square root of the non-zero eigenvalues of the covariance 
matrix 



(27) 



{Cz)kk {Cs)k,Nf + k 

iCz)k,N^+k {Cz)N^ + k,N^ + k 

zero otherwise, and t/^ is an orthogonal matrix that con- 
tains the corresponding eigenvectorqj. Indeed, if the peri- 
odogram is defined as 



Pk 



V 



Nf + kJ 



fc = 0, 1, 



1, 



(28) 



then each pk is given by the sum of two independent, unit- 
variance, Gaussian random quantities. As a consequence, 
the corresponding PDF is, independently of k, a x| whose 
cumulative distribution function (CDF) is the exponential 
function. This permits determining the statistical signifi- 
cance oipk for a specified frequency k. Things become more 
complex if Nf frequencies are inspected when looking for 
a peak. Indeed, also after the operation (pS)) . it happens 
that E[pfep;] 7^ for fc 7^ /, i.e., the frequencies of the pe- 
riodogram remain mutually correlated. This is an unavoid- 
able problem. Because of it, Nf does not correspond to 
the number of independent frequencies, so the level of false 
alarm (|20p cannot be computed. However, since Lp^ is the 
same for all the frequencies, an upper limit can be fixed for 
Lpn by setting Nf = \M/2] . The periodogram obtained by 
means of Eq. (1281) corresponds to the original Lomb-Scargle 
periodogram. 

Since the transformation ([T7| does not depend on the 
characteristics of the signal sampling, the strategy of fol- 
lowing in the case that x is the realization of (not neces- 
sarily stationary) colored noise is simply the one in Sec. [2] 
i.e. transformation of x to an array y with uncorrelated 
entries. After that, the Lomb-Scargle periodogram can be 
computed. It is worth noticing that this simple result has 
been possible thanks to a formulation of the problem in the 
time domain and the use of the matrix notation. The same 
results could have been obtained by following the popular 
approach of working in the harmonic domain but at the 
price of a much more difficult derivation. 

® The diagonal elements Ai and A2 of E* can be trivially 
computed through the solution of the quadratic equation — 
tr[C*] -|-det[C*] = 0, with tr[.] and det[.] denoting the trace and 
determinant operators. The arrays ui and 112, which constitute 
the columns of U,, can be obtained by solving the equations 
{C^-\iI)ui =0,1 = 1,2. 



4. Two examples 

To illustrate the usefulness and the simplicity of the pro- 
posed formalism in handling different situations from the 
classical ones, we show two examples in this section. 

The first consists of a periodogram of a mean-subtracted 
time series. The evaluation of the reliability of a peak in 
the periodogram of a signal x requires that (under the null 
hypothesis a; = n) n be the realization of a zero-mean 
noise process. In most experimental situations, this condi- 
tion is not fulfilled and one works with a centered (i.e. mean- 
subtracted) version x of This, however, introduces some 
(often neglected) problems. The case where x is the realiza- 
tion of a discrete white noise process with variance cr^ has 
been considered se veral times in the literature. A n exam- 
ple is the paper bv IZechmeister fc Kiirsterl (|2009[ ) where a 
rather elaborate solution is presented. With the approach 
proposed here, a simpler solution can be obtained if one 
takes into consideration that the subtraction of the mean 
from X forces a spurious correlation among the entries of x 
in such a way that the covariance matrix C-^ = E[xx^] is 
given by 

where M is number of entries of x and 1 an M x M ma- 
trix with every entry equal to unity. Since this matrix is 
singular, it cannot be diagonalized and therefore x can- 
not be whitened. In any case, if in Eq. (1241) matrix Cx is 
substituted for and one sets 



(30) 



then it is a trivial matter to decorrelate and ZN^+k by 
means of Eqs. (PS1) - (P7)) and to compute the periodogram 
through Eq. This result can be easily extend to the 

case where, because of measurement errors, each entry of x 
has its own variance cr^ . and a weighted mean is subtracted 
from the data sequence i.e. Xj — ~ J2i Vj^i/ J2i Vh with 
77/ = 1 /(T^j . Indeed, it is sufficient to substitute as given 
by Eq. ^ with 



C-^ = diag[(T^] — 



(31) 



where cr^ = [al^,crl^, 



The rest of the proce- 
dure remains unmodified. 

The second example consists of zero-mean colored noise. 
The improvement in the quality of the results obtainable 
with the approach presented in the previous section is vis- 
ible in Fig. [T] The top left panel shows a discrete signal 
Xj — 0.5 sin (27r/j) -\-nj, f = 0.127, simulated on a regu- 
lar grid of 120 time instants but with missing data in the 
ranges [31 70] and [76 115]. Here, n is the realization of 
a discrete, zero-mean, colored noise process whose auto- 
covariance function is given in the top right panel. From 
the bottom left panel, it is evident that Lomb-Scargle peri- 
odogram of the original sequence x provides rather ambigu- 
ous results concerning the presence of a sinusoidal compo- 
nent. On the other hand, such component is well visible in 
the bottom right panel that shows the Lomb-Scargle peri- 
odogram of the sequence y = C„ ' x. 
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5. Periodogram and least-squares fit of sinusoids 

The formalism proposed here is also useful in the context 
of more theoretical questions (but with important practical 
implications). For example, a point often overlooked in the 
astronomical literature is the relationship between the pe- 
riodogram and the least-squares fit of sine functions. Often 
these two methods are believed to be equivalent. Actually, 
this is true only when the sampling is regular and the fre- 
quencies of the sinusoids are given by the Fourier ones. 
Indeed, if tj = ij and fk = k/N, k = 0,1, . . . , N ~ 1, then 
Eq. (HI) can be written in the form 

X - d^a = n, (32) 

with 

and a = [ao,ai, . . . ,aN^i,bo,bi, . . . ,bN^i\'^ . The least- 
squares solution d of system (j32p is given by 

a = m'^)^dx, (34) 

wher e denotes Moore-Penrose pseudo-inverse (jBiorck I 
[199^ . In the case of even sampling, i.e. when = Fsr 
and JFx = Fx, it happens that 

d = ^x. (35) 

In other words, coefficients {uk} and {bk}, as given by 
the least-squares approach, can be obtained through the 
DFT of X, because, as shown by means of Eqs. ([TT|) . 
idd )+5 = d- In the case of uneven sampling, this identity 
is not fulfilled. Any kind of periodogram computed through 
Eq. and the least-squares fit of sine functions has to be 
expected to give different results. Moreover, as only under 
the two above-mentioned conditions do the sine functions 
constitute an orthonormal basis, the least-squares fit of a 
single sine function per time does not in general provide 
the same result as the simultaneous fit of all the sinusoids 
as in Eq. (e.g. see [Hamming |[l973l page 450). In par- 
ticular, if an unevenly-sampled signal is given by the con- 
tribution of two or more sinusoids, the one-at-a-time fit 
of a single sine function provides biased results. This also 
holds for the Lomb-Scargle periodogram, which is equiva- 
lent to the least-squares fit of a single sinusoid with a spec- 
ified frequency, with the constraint that the corresponding 
coefficients "a" and "6" are uncorrelated (Scargle . ,1983 : 
IZechmeister fc Kiirster1l2009[ ). 

6. Discussion 

As demonstrated in Sec. [3l when noise has arbitrary statis- 
tical characteristics, the computation of the periodogram 
of an unevenly-sampled signal requires two steps: 

— Whitening and standardization of the noise component 
(in this way a signal y is obtained); 

— Computation of the Lomb-Scargle periodogram of y. 

The first step, unavoidable even in the case of regular sam- 
pling, is a computationally expensive operation. Therefore, 
for time series containing more than a few thousand data 
points, dedicated algorithms exploiting the specific struc- 
ture of Cn (e.g. often this matrix is of banded type) have 



5 

to be developed for implementing Eq. (PT|). However, this 
problem is beyond the aim of the present paper. The sec- 
ond step is much less time consuming. Indeed, in the case of 
time series containing some thousands of points and when 
the periodogram has to be computed on a similar number 
of frequencies, the direct implementation of Eqs. ([^ -(^5 ]) 
results in a few seconds of computation time only. In other 
words, in many practical situations, no dedicated algorithm 
is really necessary. However, fas t algorithms have b een pro- 
posed for very long time series (|Press et al. II1992D . 

The last issue that has to be considered is to which 
extent the use of the Lomb-Scargle periodogram is really 
advantageous. Indeed, the action of the algorithms deal- 
ing with uneven sampling is essentially directed, for each 
frequency k, to force pk to be the sum of two independent 
Gaussian random quantities. However, although not clearly 
emphasized, it has already bee n pointed out t hat this op- 
eration is not critical (e.g. see IScargle I Il982l) . It can be 
expected that a periodogram computed simply through 

Pk - 2\yf, (36) 

with y = £_y, is often very close to the one given by 
Eq. (p8)) . A rigorous demonstration of this fact is difficult 
because of its strict dependence on the specific sampling. 
However, with the help of some numerical experiments and 
of the formalism introduced here, some insights are possi- 
ble. In particular, we consider the covariance matrix of 
a white noise signal when sampled on different uneven time 
grids. 

6.1. Numerical simulations 

In our simulations, we take the realization of a zero-mean, 
unit- variance, Gaussian white-noise process n sampled on 
Ms = 120 time regularly-spaced instants, but with 80 miss- 
ing data (i.e. M — 40). The available signal x can be writ- 
ten in the form 

X = Wn, (37) 
W = (Xi'Ag[w], (38) 

where W — drAg[w] and w is an array whose entries are 
equal to one in correspondence to a value of n that is avail- 
able and zero otherwise. In this case, Eq. ([3]) can be written 
as 

X = FWn (39) 
= Tn, (40) 

with = FW . Three different cases have been consid- 
ered where the missing data have time indices in the ranges 
a) [31 110], b) [31 70] and [76 115], c) [6 25], [36 75] 
and [96 115], whereas they are randomly distributed in 
a fourth case . The related covariance matrices Cj, com- 
puted through Eq. (|25p with N — Ms, are shown in Fig. [51 
whereas Fig. [3] displays the corresponding main diagonal 
of the blocks T^T^ and T^T^. Especially from Fig. [31 
it is evident that, for an arbitrary frequency k, the covari- 
ance between 5ft [Sj,] and I[£i.] is quite close to zero. This 
means that each entry of p, as given by Eq. ([55]) , can be 
assumed to be distributed according to a • From Fig. [5] 
it is also evident that, in the case of gaps present in the 
sampling pattern, there are nearby frequencies k and I for 
which not only 5ft[£j;] and 3fi[£/] are mutually correlated, but 
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also 5R[3^] with X[x_i] and 3?[x/] with These correla- 

tions, especially those between the real and the imaginary 
components, disappear in the case of random sampling. 

6.2. Interpretation of the results of the simulations 

To understand these results, it is necessary to take into 
account that Eq. p9p can be rewritten in the form 

X = FWF*Fn, (41) 
= Wn. (42) 

As is a (singular) diagonal matrix, W is a (singular) cir- 
culant matrix. This implies that x is given by the circular 
convolution of n (the DFT of the original signal, inclusive 
of the missing data) with the spectral window w (the DFT 
of the sampling pattern). Since for two arbitrary frequen- 
cies, say k and I, 2^[^fe], andl[n;] are mutually 
independent, any correlation existing between I[2fe], 
and I[x_i\ is induced by the correlation between the 
entries of di[w] and I[w]. Now, as shown in Fig. HI in the 
presence of gaps the real and the imaginary parts of w 
are both only significantly different from zero in a narrow 
interval of frequencies surrounding the origin, i.e., corre- 
sponding to the lowest frequencies {w can be interpreted 
as low-pass filter). This produces the correlations observed 
among 3fi[xj,], 3?[2/], and for nearby frequencies 
k and I. On the other hand, in the case of random sampling 
w mimics the behaviour of white noise, making the corre- 
lations less important. Moreover, Fig. [5] shows that, in the 
case of sampling with gaps, the cross-correlation between 
^[w] with I[w] is significant but not centered at zero lag. 
This explains why the quantities 3?[xj,] and k ^ I can 
be strongly correlated in spite of the small correlation be- 
tween and Similar arguments also hold for the 
case of sampling with periodic gaps (rather common in as- 
tronomical experiments). Indeed, in practical applications 
the period of these gaps is somewhat smaller than the mean 
sampling time step of an uninterrupted sequence of data. 
As a consequence, both ^[w] and I[w] present a sharp and 
narrow peak in correspondence to a frequency close to the 
origin. Actually, the spectral window of a sampling with pe- 
riodic gaps is also characterized by the presence of aliases. 
However, these aliases too are sharp and narrow, and their 
importance decreases for increasing frequencies. The com- 
bination of these facts leads again to the quantities 
3fi[x;], 3?[2;], and I[x_f.] almost being uncorrelated if the fre- 
quencies k and I are not sufficiently close enough. Moreover, 
since with periodic gaps the cross-correlation between ^[w] 
with I[w] can also be significant, but not centered on zero 
lag, then for each frequency k the correlation between 3?[xj,] 
and is negligible. These considerations are confirmed 
by Figs. [6] and [T] Figure |6] displays the covariance matrix 
Cs, computed through Eq. (^5]) . of a discrete zero-mean, 
unit- variance, and white-noise process sampled on 600 time 
instants randomly distributed (i.e. not rebinned) along six 
cycles (i.e. 100 points per cycle) of the sampling pattern 
shown in the bottom panel of Fig. [TJ 1200 frequencies have 
been considered. The top panel of Fig. [7] displays the main 
diagonal of the blocks ^sr^sr and £sr£i. 

The numerical simulations indicate that the conse- 
quences of unevenly-sampled data seem to concern the 
number of independent frequencies in p rather than the 



correlation between 3?[3j.] and its imaginary counterpart 
X[Sj,]. At present, no general method has been developed 
to deal with this problem. However, it has been pointed out 
in the literature that the number of independent frequencies 
is not a critical parameter to test the significance level of a 
peak in p. In particular, empirical arguments i ndicate that 
this number can be safely set to M/2 (e.g. see lPress et al. I 
1992). The conclusion is that forcing each entry of p to be 
the sum of two independent Gaussian quantities only has 
minor effects. This is shown by Fig. |5] where the top panel 
shows that the Lomb-Scargle periodogram of the time se- 
ries with periodic gaps considered above is quite similar to 
that provided simply by pk = 2|xj,p with 2 = ^x. This 
is evident in the bottom panel of the same figure where 
the similarity of the two periodiograms is demonstrated by 
their absolute difference. 

A final point to underline, which has important practi- 
cal implications, is that for long time series the small dif- 
ferences visible in Fig. |S] should decrease. Indeed, as seen 
above, pk will be significantly correlated with pi only if the 
two frequencies k and I are close enough. The only excep- 
tion is represented by the frequencies at the extremes of the 
frequency domain where the assumption of periodic signal 
intrinsic to DFT forces a spurious correlation. For longer 
and longer time series, this spurious correlation will affect 
a smaller and smaller fraction of frequencies and, as conse- 
quence, a larger and larger fraction of them will be mutually 
independent. This is shown in Fig. IH] where, in the context 
of the previous experiment, the mean absolute difference 
between the two periodograms is plotted as a function of 
Ns, the number of cyclic sampling patterns (Ng = 6 in 
Fig. This argument also explains why in many practi- 
cal situations the number of independent frequencies can 
be safely fixed to M/2. In conclusion, only in the case of 
signals that contain a small number of data, the Lomb- 
Scargle periodogram can be expected to exhibit noticeable 
differences from the periodogram given by Eq. Often, 
using it does not change anything. Comparable results can 
be expected with less sophisticated approaches. 

6.3. Application to an astronomical time series 

As an example of an unsophisticated method able to pro- 
duce results similar to those obtainable with the Lomb- 
Scargle periodogram, we consider the rebinning of the orig- 
inal time series on an arbitrarily dense regular time grid 
(in this way a signal with a regular sampling is obtained 
but some data are missing) followed by applying any of the 
fast Fourier transform (FFT) algorithms available nowa- 
days. Figure [10] shows an experimental (mean-subtracted) 
time series versus its rebinned version. This time series, 
which is characterized by rat her irregular samp ling, was 
obtained with the VLA array (jBiggs et al. Ill999f ) and con- 
sists of polarisation position angle measurements at an ob- 
serving frequency of 15 GHz for one of the images of the 
double gravitational lens system B0218-I-357. The original 
sequence contains only M = 45 data and it is rebinned on 
a regular grid of Mr — 92 time instants. In spite of this, 
as visible in Fig. llll the corresponding periodograms, com- 
puted on TV = Mr equispaced frequencies by means of the 
Lomb-Scargle and the FFT approach, are remarkably simi- 
lar. Here, the highest frequency approximately corresponds 
to the Nyquist frequency that is related to the shortest sam- 
pling time step. The main conclusion of this example is to 
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point out that, although in the previous section we stated 
that use of the Lomb-Scargle periodogram can be expected 
to be effective only for time series that contain a small num- 
ber of data, this is not a sufficient condition to guarantee 
that the method is truly useful. 

7. Conclusions 

In this paper wc worked out a general formalism, based on 
the matrix algebra, that is tailored to analysis of the sta- 
tistical properties of the Lomb-Scargle periodogram inde- 
pendently of the characteristics of the noise and the sam- 
pling. With this formalism it has become possible to de- 
velop a test for the presence of components of interest in 
a signal in more general situations than those considered 
in the current literature (e.g. when noise is colored and/or 
non-stationary). Moreover, we were able to clarify the rela- 
tionship between the Lomb-Scargle periodogram and other 
techniques (e.g. the least-squares fit of sinusoids) and to fix 
the conditions under which the use of such method can be 
expected to be effective. 
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Fig. 1. Results concerning the numerical experiment, presented in Sec. EJ on the detection of a sinusoidal component 
in colored noise. The top left panel shows an observed time series (blue crosses) obtained through the simulation of 
signal X = {xj}}jl^Q with Xj = 0.5 sin (27r/j) + nj, f — 0.127, on a regular grid of 120 time instants but with the 
data in the ranges [31 70] and [76 115]. Here, n (red line) is the realization of a discrete, zero-mean, colored noise 
process whose autocovariance function is given in the top right panel. For comparison, the sinusoidal component is also 
plotted (green line). The bottom left panel shows the Lomb-Scargle periodogram of the original sequence x computed 
on 120 frequencies k = 0, 1/120, 119/120, whereas the bottom right panel shows the Lomb-Scargle periodogram 
corresponding to its whitened version. In both cases, only the first 60 frequencies are shown and the vertical red line 
corresponds to the frequency of the sinusoidal component. The horizontal green line in the bottom panels provides the 
threshold corresponding to a 0.01 level of false alarm (number of independent frequencies Nf — 60), i.e. the probability 
that the periodogram of a pure noise signal exceeds such a threshold by chance is 1%. 
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Fig. 2. Covariance matrix of the real and the imaginary parts, say and of the Fourier transform {xk} of 

an unevenly sampled white-noise signal x (i.e. matrix as given by Eq. (|25p ). In particular, x is assumed to be the 
realization of a zero-mean, unit- variance, white-noise process sampled at 120 time instants regularly spaced but with 80 
missing data. Three different cases have been considered with missing data in the following ranges: a) [31 110] (top left 
panel), b) [31 70], and [76 115] (top right panel), c) [6 25], [36 75] and [96 115] (bottom left panel). In a fourth case, 
the missing data are randomly distributed (bottom right panel). Because of the small size of the figure, it is necessary 
to stress that for each panel, none of the prominent diagonal structures visible in the top right, as well in the bottom 
left quadrant, correspond to the main diagonal of the quadrant itself (i.e. none of them provide the covariance between 
3fJ[3j.] and I[x_f.] that is seen in Fig. [S]). 
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Fig. 3. Variance of and covariance of with as a function of the frequency k for the experiments in 

Fig. [5] The blue line corresponds to the main diagonal of the top left quadrant in each panel of Fig. [51 whereas the red 
one corresponds to the main diagonal in the top right, as well in the bottom left quadrant. 
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Fig. 4. Real (blue line) and imaginary (red line) parts of the spectral windows w that have been used to simulate the 
irregular sampling of the signal in the experiments corresponding to Fig. [2] (i.e., w is the Fourier transform of the sampling 
pattern w of x, see Eqs. ([571) - ([551) ). 
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Fig. 5. Cross-correlation between the real and the imaginary parts of the spectral windows w that are shown in Fig. 01 
The red cross in each panel corresponds to the point (0, 0). 
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Fig. 6. Covariance matrix of the real and the imaginary parts, say 3fJ[x;j] and of the Fourier transform {xk} of an 

unevenly sampled white- noise signal x with periodic gaps (i.e. matrix Cs as given by Eq. dH])). In particular, x is the 
realization of a discrete zero-mean, unit-variance, white-noise process sampled on a grid of 600 time instants randomly 
distributed (i.e. not rebinned) according to the cyclic sampling pattern shown in the bottom panel of Fig. [T] In total, 
1200 frequencies have been considered. 
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Fig. 7. Top panel: variance of ^[x_j.] and covariance of ^R[x_f.] with I[x_i^] (only the first 100 frequencies are shown) for the 
experiment concerning the realization of a discrete zero-mean, unit- variance, white-noise process on 600 time instants 
randomly distributed (i.e. not rebinned) according to a cyclic sampling. In total, 1200 frequencies have been considered 
(see text in Sec. Bottom panel: cyclic sampling used in the experiment. 
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Fig. 8. The top panel shows the Lomb-Scargle periodogram and the classic periodogram given by = 2|J^p for the 
experiment in Fig. [71 The bottom panel shows the smallness of their absolute difference. 
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Fig. 9. Mean absolute difference |Apfc|/A^fc (A^^ = number of frequencies) between the Lomb-Scargle periodiogram 
and the classic periodogram given by pk = 2|2j,p as a function of the number Ng of cyclic sampling patterns (Fig. [5] 
shows the case with Ng = 6). 




Fig. 10. Experimental (mean-subtracted) time series containing 45 uneven ly-spaced data versus a rebinned version 
computed on a regular grid of 92 time instants (data taken from .Biggs et al. Ill999i) . 
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Fig. 11. Periodograms of the time series shown in Fig. 1101 Frequency is in given in units of the Nyquist frequency 
corresponding to the shortest samphng time step. The periodogram of the original time series has been obtained by 
means of the Lomb-Scargle method and that of the rebinned version by means of a classic FFT algorithm. 



